CSE 252A Computer Vision I Fall 2019 - Homework 4¶
Instructor: Ben Ochoa¶
Assignment published on: Tuesday, November 05, 2019¶
Due on: Tuesday, November 19, 2019 11:59 pm¶
Instructions¶
- Review the academic integrity and collaboration policies on the course website.
- This assignment must be completed individually.
- All solutions must be written in this notebook.
- This includes the theoretical problems, for which you must write your answers in Markdown cells (using LaTeX when appropriate).
- Programming aspects of the assignment must be completed using Python in this notebook.
- If you want to modify the skeleton code, you may do so. It has only been provided as a framework for your solution.
- You may use Python packages (such as NumPy and SciPy) for basic linear algebra, but you may not use packages that directly solve the problem.
- If you are unsure about using a specific package or function, then ask the instructor and/or teaching assistants for clarification.
- You must submit this notebook exported as a PDF. You must also submit this notebook as
.ipynbfile.- Submit both files (
.pdfand.ipynb) on Gradescope. - You must mark the PDF pages associated with each question in Gradescope. If you fail to do so, we may dock points.
- Submit both files (
- It is highly recommended that you begin working on this assignment early.
- Late policy: assignments submitted late will receive a 15% grade reduction for each 12 hours late (i.e., 30% per day). Assignments will not be accepted 72 hours after the due date. If you require an extension (for personal reasons only) to a due date, you must request one as far in advance as possible. Extensions requested close to or after the due date will only be granted for clear emergencies or clearly unforeseeable circumstances.
Problem 1: Epipolar Geometry [4 pts]¶
Consider two cameras whose image planes are the z=1 plane, and whose focal points are at (-12, 0, 0) and (12, 0, 0). We''ll call a point in the first camera (x, y), and a point in the second camera (u, v). Points in each camera are relative to the camera center. So, for example if (x, y) = (0, 0), this is really the point (-12, 0, 1) in world coordinates, while if (u, v) = (0, 0) this is the point (12, 0, 1).
a) Suppose the points (x, y) = (8, 7) is matched to the point (u, v) = (2, 7). What is the 3D location of this point?
Solution: Given that, Baseline $d = 12-(-12) = 24$, the positions of camera1 and camera2 are corespondingly $C_1 = (-\frac{d}{2},0,0)$ and $C_2 = (\frac{d}{2},0,0)$, then use similar triangle, we have, $$ \begin{cases} \frac{X_L}{X-(-d/2)} = \frac{f}{Z}\\ \frac{X_R}{X-{d/2}} = \frac{f}{Z}\\ \frac{Y}{Z} = \frac{y}{f} \end{cases} $$ solving the equations, we get$\begin{cases} X = 20\\ Y = 28 \\Z = 4 \end{cases}$, thus the 3D location of this point is (20,28,4).
b) Compute the Essential Matrix.
Solution:
Since the image plane has not been rotated, the rotation matrix $R$ should be an identical matrix, and the transform matrix $t$ from Camera1 to Camera2 is along the x-axis, so $$ R = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ $$ t = [t_x, t_y, t_z]^T \Rightarrow [t_{\times}] = \begin{bmatrix} 0 & -t_z & t_y \\ t_z & 0 & -t_x \\ -t_y & t_x & 0 \end{bmatrix} $$ Thus, the essential matrix $$ E = R[t_{\times}] = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -24 \\ 0 & 24 & 0 \end {bmatrix} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -24 \\ 0 & 24 & 0 \end{bmatrix} $$
c) Consider points that lie on the line x + z = 2, y = 0. Use the same stereo set up as before. Write an analytic expression giving the disparity of a point on this line after it projects onto the two images, as a function of its position in the right image. So your expression should only involve the variables u and d (for disparity). Your expression only needs to be valid for points on the line that are in front of the cameras, i.e. with z > 1.
Solution:
by definition, disparity $d = X_L - X_R$
By similar triangle, we have $$ \frac{d}{24}=\frac{f}{z}=\frac{u}{12-x} $$ with a constrain of $$ x+z=2 $$ Thus, $$ d = 2.4u-2.4 $$
Problem 2: Epipolar Rectification [4 pts]¶
In stereo vision, image rectification is a common preprocessing step to simplify the problem of finding matching points between images. The goal is to warp image views such that the epipolar lines are horizontal scan lines of the input images. Suppose that we have captured two images $I_A$ and $I_B$ from identical calibrated cameras separated by a rigid transformation
$_{A}^{B}\textrm{T}= \begin{bmatrix} \boldsymbol{R} & \boldsymbol{t} \\ 0^T & 1 \end{bmatrix}$ and $\boldsymbol{t}=[t_x,t_y,t_z]$
Without loss of generality assume that camera A's optical center is positioned at the origin and that its optical axis is in the direction of the z-axis.
From the lecture, a rectifying transform for each image should map the epipole to a point infinitely far away in the horizontal direction $ H_{A}e_{A} = H_{B}e_{B} = [1, 0, 0]^T$. Consider the following special cases:
a) Pure horizontal translation $\boldsymbol{t} = [t_{x}, 0, 0]^T$, $\boldsymbol{R} = \boldsymbol{I}$
b) Pure translation orthogonal to the optical axis $\boldsymbol{t} = [t_{x}, t_{y}, 0]^T$, $\boldsymbol{R} = \boldsymbol{I}$
c) Pure translation along the optical axis $\boldsymbol{t} = [0, 0, t_{z}]^T$, $\boldsymbol{R} = \boldsymbol{I}$
d) Pure rotation $\boldsymbol{t} = [0, 0, 0]^T$, $\boldsymbol{R}$ is an arbitrary rotation matrix
For each of these cases, determine whether or not epipolar rectification is possible. Include the following information for each case:
(i) The epipoles $e_A$ and $e_B$
(ii) The equation of the epipolar line $l_B$ in $I_B$ corresponding to the point $[x_A, y_A, 1]^T$ in $I_A$ (if one exists)
(iii) A plausible solution to the rectifying transforms $H_A$ and $H_B$ (if one exists) that attempts to minimize distortion (is as close as possible to a 2D rigid transformation). Note that the above 4 cases are special cases; a simple solution should become apparent by looking at the epipolar lines.
One or more of the above rigid transformations may be a degenerate case where rectification is not possible or epipolar geometry does not apply. If so, explain why.
Solution:
a) Given that $R = I$, ${t} = [t_{x}, 0, 0]^T$, $$ E = R[t_{\times}] = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -t_x \\ 0 & t_x & 0 \end{bmatrix} $$ i) $$ \begin{cases} E\cdot e_B = 0\\ E^T\cdot e_A = 0 \end{cases} \Rightarrow e_A=e_B = [1, 0, 0]^T $$ ii) $$ [a,b,c]^T = E^T\cdot[x_A,y_A,1]^T=[0,t_x,-t_x y_A]^T $$ Thus, $l_B$ is $$ t_x\cdot y - t_x \cdot y_A = 0 \Rightarrow y=y_A $$ iii) $$ e_1 = e_A = [1,0,0], e_2 = \frac{1}{\sqrt{e_{1x}^2+e_{1y}^2}}[-e_{1y},e_{1x},0]^T=[0,1,0]^T, e_3 = e_1 \times e_2 = [0,0,1]^T $$ Thus, $$ H_A = [e_1^T, e_2^T, e_3^T]^T = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, H_B = R\cdot H_A = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ To conclude, epipolar rectification is possible.
b) For $R = I$,${t} = [t_{x}, t_{y}, 0]^T$, we have, $$ E = R[t_{\times}] = \begin{bmatrix} 0 & 0 & t_y \\ 0 & 0 & -t_x \\ -t_y & t_x & 0 \end{bmatrix} $$ i) $$ \begin{cases} E\cdot e_B = 0\\ E\cdot e_A = 0 \end{cases} \Rightarrow e_A = e_B = \left[\frac{t_x}{\sqrt{t_x^2+t_y^2}},\frac{t_y}{\sqrt{t_x^2+t_y^2}},0\right]^T $$ ii) $$ [a,b,c] = E^T[x_A,y_A,1]^T=[-t_y, t_x,t_y x_A-t_x y_A]^T $$ Thus, $l_B$ is $$ -t_y\cdot x + t_x\cdot y + t_y x_A - t_x y_A = 0 $$ iii) $$ e_1 = e_A,\\ e_2 = \frac{1}{\sqrt{e_{1x}^2+e_{1y}^2}}[-e_{1y},e_{1x},0]^T = \left[ -\frac{t_y}{\sqrt{t_x^2+t_y^2}},\frac{t_x}{\sqrt{t_x^2+t_y^2}},0 \right]^T\\ e_3 = e_1 \times e_2 = [0,0,1]^T $$ Thus, $$ H_A = \begin{bmatrix} e_1^T\\ e_2^T\\ e_3^T \end{bmatrix} =\begin{bmatrix} \frac{t_x}{\sqrt{t_x^2+t_y^2}} & \frac{t_y}{\sqrt{t_x^2+t_y^2}} & 0 \\ -\frac{t_y}{\sqrt{t_x^2+t_y^2}} & \frac{t_x}{\sqrt{t_x^2+t_y^2}} & 0\\ 0 & 0 & 0 \end{bmatrix} =R\cdot H_A=H_B $$ To conclude, epipolar rectification is possible.
c) Given that $R = I$, ${t} = [0, 0,t_{y}]^T$, we have, $$ E = R[t_{\times}] = \begin{bmatrix} 0 & -t_z & 0 \\ t_z & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$ i) $$ \begin{cases} E\cdot e_B = 0\\ E\cdot e_A = 0 \end{cases} \Rightarrow e_A = e_B =[0,0,1]^T $$ ii) $$ [a,b,c]^T = E^T\cdot [x_A,y_A,1]^T=[t_z y_A, -t_z x_A, 0]^T $$ Thus, $l_B$ can be written as: $$ t_z y_A \cdot x - t_z x_A\cdot y = 0 \Rightarrow y_A \cdot x - x_A\cdot y=0 $$ iii) Since $e_1=e_A=[0,0,1]^T$, suppose $e_2=[0,1,0]^T$ s.t. $e_2 \perp e_1$.
Thus, $e_3=e_1\times e_2 = [1,0,0]^T$, and then we have $$ H_A = \begin{bmatrix} 0&0&1 \\ 0&1&0 \\ 1&0&0 \end{bmatrix}\quad and \quad H_B = R\cdot H_A = \begin{bmatrix} 0&0&1 \\ 0&1&0 \\ 1&0&0 \end{bmatrix} $$ To conclude, under this situation, we can randomly pick a unit vector $e_2$ on x-y plane such that this vector $e_2$ is orthogonal to vector $e_1$ and then generate homography matrix. Epipolar rectification is possible.
d) Given that ${t} = [0,0,0]^T$, $\forall R$, we have $$ E = R[t_{\times}] = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$ Obviously, we cannot pick a specific unit vector $e_1$ since the rotation is arbitrary. As a result of this, $l_B$ does not exist, so does $H_A$ and $ H_B$.
Thus, epipolar rectification is impossible.
Problem 3: Sparse Stereo Matching [32 pts]¶
In this problem we will play around with sparse stereo matching methods. You will work on two image pairs, a warrior figure and a figure from the Matrix movies. These files both contain two images, two camera matrices, and two sets of corresponding points (extracted by manually clicking the images). For illustration, I have run my code on a third image pair (dino1.png, dino2.png). This data is also provided for you to debug your code, but you should only report results on warrior and matrix. In other words, where I include one (or a pair) of images in the assignment below, you will provide the same thing but for BOTH matrix and warrior. Note that the matrix image pair is harder, in the sense that the matching algorithms we are implementing will not work quite as well. You should expect good results, however, on warrior.
Corner Detection [8 pts]¶
The first thing we need to do is to build a corner detector. This should be done according to http://cseweb.ucsd.edu/classes/fa19/cse252A-a/lec7.pdf. You should fill in the function corner_detect below, and take as input corner_detect(image, nCorners, smoothSTD, windowSize) where smoothSTD is the standard deviation of the smoothing kernel and windowSize is the window size for corner detector and non maximum suppression. In the lecture the corner detector was implemented using a hard threshold. Do not do that but instead return the nCorners strongest corners after non-maximum suppression. This way you can control exactly how many corners are returned. Run your code on all four images (with nCorners = 20) and show outputs as shown below. You may find scipy.ndimage.filters.gaussian_filter easy to use for smoothing. In this problem, try the following different standard deviation ($\sigma$) parameters for the Gausian smoothing kernel: 0.5, 1, 2 and 4. For a particular $\sigma$, you should take the kernel size to be $6\times\sigma$ (add $1$ if the kernel size is even). So for example if $\sigma=2$, corner detection kernel size should be $13$. This should be followed throughtout all experiments in this assignment.
There will be a total of 16 images as outputs : (4 choices of smoothSTD x 2 matrix imgs x 2 warrior imgs).
Comment on your results and observations (3/8 points). You don't need to comment per output, just discuss any trends you see for the detected corners as you change the windowSize and increase the smoothing w.r.t the two pair of images warrior and matrix. Also discuss if you are able to find corresponding corners for same pair of images.
import numpy as np
from scipy.misc import imread
from scipy.signal import convolve
import matplotlib.pyplot as plt
from scipy.ndimage.filters import gaussian_filter
import imageio
def rgb2gray(rgb):
""" Convert rgb image to grayscale.
"""
return np.dot(rgb[...,:3], [0.299, 0.587, 0.114])
def corner_detect(image, nCorners, smoothSTD, windowSize):
"""Detect corners on a given image.
Args:
image: Given a grayscale image on which to detect corners.
nCorners: Total number of corners to be extracted.
smoothSTD: Standard deviation of the Gaussian smoothing kernel.
windowSize: Window size for corner detector and non maximum suppression.
Returns:
Detected corners (in image coordinate) in a numpy array (n*2).
"""
"""
Put your awesome numpy powered code here:
"""
radi = windowSize // 2
img_smth = gaussian_filter(image, sigma=smoothSTD)
dx_kernel = np.array([[-0.5, 0, 0.5]])
dx_img = convolve(img_smth, dx_kernel, mode='same')
dx_img[:, 0] = dx_img[:, 1]
dx_img[:, -1] = dx_img[:, -2]
dy_kernel = np.array([[-0.5, 0, 0.5]]).T
dy_img = convolve(img_smth, dy_kernel, mode='same')
dy_img[0, :] = dy_img[1, :]
dy_img[-1, :] = dy_img[-2, :]
C_lambda = np.zeros([image.shape[0], image.shape[1]])
dxx = dx_img*dx_img
dxy = dx_img*dy_img
dyy = dy_img*dy_img
for row in range(image.shape[0]):
for col in range(image.shape[1]):
top = 0 if (row - radi < 0) else row - radi
bottom = image.shape[0] if (radi + row > image.shape[0]) else radi + row
left = 0 if (col - radi) < 0 else col - radi
right = image.shape[1] if (radi + col > image.shape[1]) else radi + col
Ix2 = np.sum(dxx[top:bottom + 1, left:right + 1])
Ixy = np.sum(dxy[top:bottom + 1, left:right + 1])
Iy2 = np.sum(dyy[top:bottom + 1, left:right + 1])
c = np.array([[Ix2, Ixy], [Ixy, Iy2]])
C_lambda[row, col] = min(np.linalg.eigvals(c))
# nms
C_nms = np.array([0, 0, 0])
for row in range(image.shape[0]):
for col in range(image.shape[1]):
top = 0 if (row - radi < 0) else row - radi
bottom = image.shape[0] if (radi + row > image.shape[0]) else radi + row
left = 0 if (col - radi) < 0 else col - radi
right = image.shape[1] if (radi + col > image.shape[1]) else radi + col
cWindow = C_lambda[top:bottom + 1, left:right + 1]
maxLambda = max(cWindow.flatten())
if(maxLambda == C_lambda[row,col]):
C_nms = np.vstack((C_nms,np.array([maxLambda, col, row])))
else:
continue
C_nms = np.unique(C_nms, axis=0)
C_nms_sort = C_nms[np.lexsort(-C_nms[:, ::-1].T)]
corners = np.zeros((nCorners, 2))
for rowCorner in range(nCorners):
corners[rowCorner][0] = C_nms_sort[rowCorner][1]
corners[rowCorner][1] = C_nms_sort[rowCorner][2]
return corners
def show_corners_result(imgs, corners):
fig = plt.figure(figsize=(15, 15))
ax1 = fig.add_subplot(221)
ax1.imshow(imgs[0], cmap='gray')
ax1.scatter(corners[0][:, 0], corners[0][:, 1], s=36, edgecolors='r', facecolors='none')
ax2 = fig.add_subplot(222)
ax2.imshow(imgs[1], cmap='gray')
ax2.scatter(corners[1][:, 0], corners[1][:, 1], s=36, edgecolors='r', facecolors='none')
plt.show()
# detect corners on warrior and matrix image sets
# adjust your corner detection parameters here
nCorners = 20
smoothSTDs = [0.5, 1, 2, 4]
imgs_mat = []
imgs_war = []
grayimgs_mat = []
grayimgs_war = []
# Read the two images and convert it to Greyscale
for i in range(2):
img_mat = imageio.imread('p4/matrix/matrix' + str(i) + '.png')
imgs_mat.append(img_mat)
grayimgs_mat.append(rgb2gray(img_mat))
# Comment above line and uncomment below line to
# downsize your image in case corner_detect runs slow in test
# grayimgs_mat.append(rgb2gray(img_mat)[::2, ::2])
# if you unleash the power of numpy you wouldn't need to downsize, it'll be fast
img_war = imageio.imread('p4/warrior/warrior' + str(i) + '.png')
imgs_war.append(img_war)
grayimgs_war.append(rgb2gray(img_war))
# grayimgs_war.append(rgb2gray(img_war)[::2,::2])
for smoothSTD in smoothSTDs:
windowSize = int(6*smoothSTD)
if windowSize%2==0: windowSize += 1
crns_mat = []
crns_war = []
print ("SmoothSTD:", smoothSTD, "WindowSize:", windowSize)
for i in range(2):
crns_mat.append(corner_detect(grayimgs_mat[i], nCorners, smoothSTD,\
windowSize))
crns_war.append(corner_detect(grayimgs_war[i], nCorners, smoothSTD,\
windowSize))
show_corners_result(imgs_mat, crns_mat) #uncomment this to show your output!
show_corners_result(imgs_war, crns_war)
With the increasing of SmoothSTD
Because of the texture of clothes of the figure in matrix.png, the figure's chest part and the sleeves part can be misidentified as corner when the contrast is obvious. The gaussian kernel can decrease the contrast of each adjunct pixel, after using a larger scale of gaussian smoothing kernel, the image becomes "smoother".
It can be observed that with the increasing of the value of smoothSTD, the corners are more likely to be detected along the outskirt of the figure and show up less frequently on the clothes of that figure. And this implement makes the corner detection more accurate.
With the increasing of WindowSize
The WindowSize are utilized in calculating the eigenvalues of the windows and implementing the NMS. The corners of the warrior figure can be more concentrate on the forehead when the WindowSize is small. Since a larger WindowSize guarantees a detection of a "larger" corner, the eigenvalue of the window around the figure's forehead area will no longer be higher than other actual corners.
Also, using the NMS method, the corners which are too close to each other will be suppressed. So, the corners will be more widely distributed on around the figure. The larger the WindowSize is, the more seperate distributed the corners are.
Can the corners be matched?
Under the same enviornment light source, the paired pictures should look similar in general. If the corners from both images, they should have similar "calculate features". For example, similar eigenvalues in the corresponding window, or the highest "matching score" that we implemented. Those are some features for distinguishing whether the corresponding corners match each other or not.
NCC (Normalized Cross-Correlation) Matching [2 pts]¶
Write a function ncc_match that implements the NCC matching algorithm for two input windows.
NCC = $\sum_{i,j}\tilde{W_1} (i,j)\cdot \tilde{W_2} (i,j)$ where $\tilde{W} = \frac{W - \overline{W}}{\sqrt{\sum_{k,l}(W(k,l) - \overline{W})^2}}$ is a mean-shifted and normalized version of the window and $\overline{W}$ is the mean pixel value in the window W.
def ncc_match(img1, img2, c1, c2, R):
"""Compute NCC given two windows.
Args:
img1: Image 1.
img2: Image 2.
c1: Center (in image coordinate) of the window in image 1.
c2: Center (in image coordinate) of the window in image 2.
R: R is the radius of the patch, 2 * R + 1 is the window size
Returns:
NCC matching score for two input windows.
"""
"""
Your code here:
"""
matching_score = 0
w1_left,w1_right = c1[1]-R, c1[1]+R+1
w1_top, w1_bottom = c1[0]-R, c1[0]+R+1
w2_left,w2_right = c2[1]-R, c2[1]+R+1
w2_top, w2_bottom = c2[0]-R, c2[0]+R+1
window1 = img1[w1_left:w1_right,w1_top:w1_bottom]
window2 = img2[w2_left:w2_right,w2_top:w2_bottom]
W1_mean = np.mean(window1)
W2_mean = np.mean(window2)
temp1 = np.sqrt(np.sum(np.square(window1-W1_mean)))
temp2 = np.sqrt(np.sum(np.square(window2-W2_mean)))
for row in range(window1.shape[0]):
for col in range(window1.shape[1]):
w1_temp = (window1[row,col]-W1_mean)/temp1
w2_temp = (window2[row,col]-W2_mean)/temp2
matching_score += w1_temp*w2_temp
return matching_score
# test NCC match
img1 = np.array([[1, 2, 3, 4], [4, 5, 6, 8], [7, 8, 9, 4]])
img2 = np.array([[1, 2, 1, 3], [6, 5, 4, 4], [9, 8, 7, 3]])
print (ncc_match(img1, img2, np.array([1, 1]), np.array([1, 1]), 1))
# should print 0.8546
print (ncc_match(img1, img2, np.array([2, 1]), np.array([2, 1]), 1))
# should print 0.8457
print (ncc_match(img1, img2, np.array([1, 1]), np.array([2, 1]), 1))
# should print 0.6258
Naive Matching [4 pts]¶
Equipped with the corner detector and the NCC matching function, we are ready to start finding correspondances. One naive strategy is to try and find the best match between the two sets of corner points. Write a script that does this, namely, for each corner in image1, find the best match from the detected corners in image2 (or, if the NCC match score is too low, then return no match for that point). You will have to figure out a good threshold (NCCth) value by experimentation. Write a function naive_matching and call it as below. Examine your results for 10, 20, and 30 detected corners in each image. Choose number of detected corners to maximize the number of correct matching pairs. naive_matching will call your NCC matching code.
Properly label or mention which output corresponds to which choice of number of corners. Total number of output is 6 images (3 choice of number of corners for each matrix and warrior), where one image is like below:
Number of Corners: 10
def pointsMatch(matrix):
"""
***Feature matching***
calculate correlation coefficient between a
given window in image1 and all windows in image2
(window about corner point)
correlation coeffient value is [-1,1] (Greater is better)
***One-to-one matching***
1. Find indices of the element with maximum value in masked correlation
coefficient array, this is the best match.(now check if it is unique enough)
2. store best match value
3. temporarily, set the value of the element in
the correlation coefficient array to -1
4. Find next best match value as
Max(max value in same row as element,
max value in same column of element) not masked
5. set the value of the element back to its oringinal value
6. if (1 - best match value)<(1 - next best match value)* distance ratio threshold:
store feature match
(else, match is not unique enough)
7. in mask array, set row and column corresponding to best match to false
8. repeat until similarity threshold >= maximum value in masked correlation
coefficient array(determined in step 1)
"""
if np.amax(matrix) == 0:
return -matrix
arg = np.argmax(matrix)
max_tp = np.amax(matrix)
matrix[arg//matrix.shape[0],arg-(arg//matrix.shape[0]+1)*matrix.shape[0]]= -1
row_tp = matrix[arg//matrix.shape[0],:]
col_tp = matrix[:,arg-(arg//matrix.shape[0]+1)*matrix.shape[0]]
max2nd = max(np.amax(row_tp),np.amax(col_tp))
if (1-max_tp)<((1-max2nd)*0.8):
for i in range(matrix.shape[0]):
matrix[arg//matrix.shape[0],i] = 0
matrix[i,arg-(arg//matrix.shape[0]+1)*matrix.shape[0]] = 0
matrix[arg//matrix.shape[0],arg-(arg//matrix.shape[0]+1)*matrix.shape[0]]= -1
else:
for i in range(matrix.shape[0]):
matrix[arg//matrix.shape[0],i] = 0
matrix[i,arg-(arg//matrix.shape[0]+1)*matrix.shape[0]] = 0
return pointsMatch(matrix)
def naive_matching(img1, img2, corners1, corners2, R, NCCth):
"""Compute NCC given two windows.
Args:
img1: Image 1.
img2: Image 2.
corners1: Corners in image 1 (nx2)
corners2: Corners in image 2 (nx2)
R: NCC matching radius
NCCth: NCC matching score threshold
Returns:
NCC matching result a list of tuple (c1, c2),
c1 is the 1x2 corner location in image 1,
c2 is the 1x2 corner location in image 2.
"""
"""
Your code here:
"""
matching = []
score_mat = np.zeros([corners1.shape[0],corners2.shape[0]])
ncc_max = 0
for c1 in range(corners1.shape[0]):
for c2 in range(corners2.shape[0]):
score_12 = ncc_match(img1, img2, np.array(corners1[c1],dtype=int),np.array(corners2[c2],dtype=int),R)
if score_12 < NCCth :
continue
score_mat[c1,c2] = score_12
score_mat = pointsMatch(score_mat)
for row in range(score_mat.shape[0]):
for col in range(score_mat.shape[1]):
if score_mat[row,col] == 0:
continue
else:
matching.append((corners1[row,:],corners2[col,:]))
return matching
# You are free to modify code here, create your helper functions etc.
# detect corners on warrior and matrix sets
nCorners = 30 # Do this for 10, 20 and 30 corners
smoothSTD = 2
windowSize = 13
# read images and detect corners on images
imgs_mat = []
crns_mat = []
imgs_war = []
crns_war = []
for i in range(2):
img_mat = imageio.imread('p4/matrix/matrix' + str(i) + '.png')
imgs_mat.append(rgb2gray(img_mat))
# downsize your image in case corner_detect runs slow in test
# imgs_mat.append(rgb2gray(img_mat)[::2, ::2])
crns_mat.append(corner_detect(imgs_mat[i], nCorners, smoothSTD, windowSize))
img_war = imageio.imread('p4/warrior/warrior' + str(i) + '.png')
imgs_war.append(rgb2gray(img_war))
# imgs_war.append(rgb2gray(img_war)[::2, ::2])
crns_war.append(corner_detect(imgs_war[i], nCorners, smoothSTD, windowSize))
# match corners
crns_mat_10_1 = crns_mat[0][0:10]
crns_mat_10_2 = crns_mat[1][0:10]
crns_mat_10 = []
crns_mat_10.append(crns_mat_10_1)
crns_mat_10.append(crns_mat_10_2)
crns_mat_20_1 = crns_mat[0][0:20]
crns_mat_20_2 = crns_mat[1][0:20]
crns_mat_20 = []
crns_mat_20.append(crns_mat_20_1)
crns_mat_20.append(crns_mat_20_2)
crns_war_10_1 = crns_war[0][0:10]
crns_war_10_2 = crns_war[1][0:10]
crns_war_10 = []
crns_war_10.append(crns_war_10_1)
crns_war_10.append(crns_war_10_2)
crns_war_20_1 = crns_war[0][0:20]
crns_war_20_2 = crns_war[1][0:20]
crns_war_20 = []
crns_war_20.append(crns_war_20_1)
crns_war_20.append(crns_war_20_2)
R = 15
NCCth = 0.7 # Put your threshold
matching_mat10 = naive_matching(imgs_mat[0]/255, imgs_mat[1]/255, crns_mat_10[0], crns_mat_10[1], R, NCCth)
matching_war10 = naive_matching(imgs_war[0]/255, imgs_war[1]/255, crns_war_10[0], crns_war_10[1], R, NCCth)
matching_mat20 = naive_matching(imgs_mat[0]/255, imgs_mat[1]/255, crns_mat_20[0], crns_mat_20[1], R, NCCth)
matching_war20 = naive_matching(imgs_war[0]/255, imgs_war[1]/255, crns_war_20[0], crns_war_20[1], R, NCCth)
matching_mat30 = naive_matching(imgs_mat[0]/255, imgs_mat[1]/255, crns_mat[0], crns_mat[1], R, NCCth)
matching_war30 = naive_matching(imgs_war[0]/255, imgs_war[1]/255, crns_war[0], crns_war[1], R, NCCth)
# plot matching result
def show_matching_result(img1, img2, matching):
fig = plt.figure(figsize=(8, 8))
plt.imshow(np.hstack((img1, img2)), cmap='gray') # two dino images are of different sizes, resize one before use
for p1, p2 in matching:
plt.scatter(p1[0], p1[1], s=35, edgecolors='r', facecolors='none')
plt.scatter(p2[0] + img1.shape[1], p2[1], s=35, edgecolors='r', facecolors='none')
plt.plot([p1[0], p2[0] + img1.shape[1]], [p1[1], p2[1]])
plt.savefig('dino_matching.png')
plt.show()
# Uncomment to show output
print("Number of Corners:", 10)
show_matching_result(imgs_mat[0], imgs_mat[1], matching_mat10)
show_matching_result(imgs_war[0], imgs_war[1], matching_war10)
print("Number of Corners:", 20)
show_matching_result(imgs_mat[0], imgs_mat[1], matching_mat20)
show_matching_result(imgs_war[0], imgs_war[1], matching_war20)
print("Number of Corners:", 30)
show_matching_result(imgs_mat[0], imgs_mat[1], matching_mat30)
show_matching_result(imgs_war[0], imgs_war[1], matching_war30)
Epipolar Geometry [4 pts]¶
Complete the compute_fundamental function below using 8 point algorithm described in Lecture 8. Using the fundamental_matrix function and the corresponding points provided in cor1.npy and cor2.npy, calculate the fundamental matrix for the set of matrix and warrior image. Note that the normalization of the corner point is handled in the fundamental_matrix function.
import numpy as np
from scipy.misc import imread
import matplotlib.pyplot as plt
from scipy.io import loadmat
def compute_fundamental(x1,x2):
""" Computes the fundamental matrix from corresponding points
(x1,x2 3*n arrays) using the 8 point algorithm.
Each row in the A matrix below is constructed as
[x'*x, x'*y, x', y'*x, y'*y, y', x, y, 1]
Returns:
Fundamental Matrix (3x3)
"""
"""
Your code here
"""
n = x1.shape[1]
if x2.shape[1] != n:
raise ValueError("Number of points don't match.")
A = np.zeros((n,9))
for i in range(n):
A[i] = [x1[0,i]*x2[0,i], x1[0,i]*x2[1,i], x1[0,i],
x1[1,i]*x2[0,i], x1[1,i]*x2[1,i], x1[1,i],
x2[0,i], x2[1,i], 1 ]
u, sigma, v = np.linalg.svd(A)
f = v[-1].reshape((3,3),order='C')
uf, sigmaf,vf = np.linalg.svd(f)
sigmaf[2] = 0
F = np.dot(uf,np.dot(np.diag(sigmaf),vf))
return F/F[2,2]
def fundamental_matrix(x1,x2):
# Normalization of the corner points is handled here
n = x1.shape[1]
if x2.shape[1] != n:
raise ValueError("Number of points don't match.")
# normalize image coordinates
x1 = x1 / x1[2]
mean_1 = np.mean(x1[:2],axis=1)
S1 = np.sqrt(2) / np.std(x1[:2])
T1 = np.array([[S1,0,-S1*mean_1[0]],[0,S1,-S1*mean_1[1]],[0,0,1]])
x1 = np.dot(T1,x1)
x2 = x2 / x2[2]
mean_2 = np.mean(x2[:2],axis=1)
S2 = np.sqrt(2) / np.std(x2[:2])
T2 = np.array([[S2,0,-S2*mean_2[0]],[0,S2,-S2*mean_2[1]],[0,0,1]])
x2 = np.dot(T2,x2)
# compute F with the normalized coordinates
F = compute_fundamental(x1,x2)
# reverse normalization
F = np.dot(T1.T,np.dot(F,T2))
return F/F[2,2]
Plot Epipolar Lines [5 pts]¶
Using this fundamental matrix, plot the epipolar lines in both image pairs across all images. For this part you may want to complete the function plot_epipolar_lines. Shown your result for matrix and warrior as the figure below.
Also, write the script to calculate the epipoles for a given Fundamental matrix and corner point correspondences in the two images.
def compute_epipole(F):
'''
This function computes the epipoles for a given fundamental matrix
and corner point correspondences
input:
F--> Fundamental matrix
output:
e1--> corresponding epipole in image 1
e2--> epipole in image2
'''
#your code here
u1,s1,v1 = np.linalg.svd(F.T)
e1 = v1[-1,:]/v1[-1,:][-1]
u2,s2,v2 = np.linalg.svd(F)
e2 = v2[-1,:]/v2[-1,:][-1]
return e1,e2
def plot_epipolar_lines(img1,img2, cor1, cor2):
"""Plot epipolar lines on image given image, corners
Args:
img1: Image 1.
img2: Image 2.
cor1: Corners in homogeneous image coordinate in image 1 (3xn)
cor2: Corners in homogeneous image coordinate in image 2 (3xn)
"""
"""
Your code here:
"""
F = fundamental_matrix(cor1, cor2)
e1,e2 = compute_epipole(F)
fig2 = plt.figure(figsize=(8, 8))
plt.imshow(img1)
for i in range(cor1.shape[1]):
plt.scatter(cor1[0][i],cor1[1][i],s=35, edgecolors='b',facecolors='b')
for i in range(cor1.shape[1]):
k_cor1=e1-cor1[:,i]
k_cor1=k_cor1[1]/k_cor1[0]
y1 = k_cor1*(0-cor1[0,i])+cor1[1,i]
y2 = k_cor1*(img1.shape[1]-cor1[0,i])+cor1[1,i]
plt.plot([0, img1.shape[1]], [y1, y2], color = 'b')
plt.imshow(img1)
plt.show()
fig2 = plt.figure(figsize=(8, 8))
plt.imshow(img2)
for i in range(cor2.shape[1]):
plt.scatter(cor2[0][i],cor2[1][i],s=35, edgecolors='b',facecolors='b')
for i in range(cor2.shape[1]):
k_cor2=e2-cor2[:,i]
k_cor2=k_cor2[1]/k_cor2[0]
y1 = k_cor2*(0-cor2[0,i])+cor2[1,i]
y2 = k_cor2*(img1.shape[1]-cor2[0,i])+cor2[1,i]
plt.plot([0,img2.shape[1]], [y1, y2], color = 'b')
plt.imshow(img2)
plt.show()
# replace images and corners with those of matrix and warrior
imgids = ["dino", "matrix", "warrior"]
for imgid in imgids:
I1 = imageio.imread("./p4/"+imgid+"/"+imgid+"0.png")
I2 = imageio.imread("./p4/"+imgid+"/"+imgid+"1.png")
cor1 = np.load("./p4/"+imgid+"/cor1.npy")
cor2 = np.load("./p4/"+imgid+"/cor2.npy")
plot_epipolar_lines(I1,I2,cor1,cor2)
Image Rectification [5 pts]¶
An interesting case for epipolar geometry occurs when two images are parallel to each other. In this case, there is no rotation component involved between the two images and the essential matrix is $\texttt{E}=[\boldsymbol{T_{x}}]\boldsymbol{R}=[\boldsymbol{T_{x}}]$. Also if you observe the epipolar lines $\boldsymbol{l}$ and $\boldsymbol{l^{'}}$ for parallel images, they are horizontal and consequently, the corresponding epipolar lines share the same vertical coordinate. Therefore the process of making images parallel becomes useful while discerning the relationships between corresponding points in images.
Rectifying a pair of images can also be done for uncalibrated camera images (i.e. we do not require the camera matrix of intrinsic parameters). Using the fundamental matrix we can find the pair of epipolar lines $\boldsymbol{l_i}$ and $\boldsymbol{l^{'}_i}$ for each of the correspondences. The intersection of these lines will give us the respective epipoles $\boldsymbol{e}$ and $\boldsymbol{e^{'}}$. Now to make the epipolar lines to be parallel we need to map the epipoles to infinity. Hence , we need to find a homography that maps the epipoles to infinity. The method to find the homography has been implemented for you. You can read more about the method used to estimate the homography in the paper "Theory and Practice of Projective Rectification" by Richard Hartley.
Using the compute_epipoles function from the previous part and the given compute_matching_homographies function, find the rectified images and plot the parallel epipolar lines using the plot_epipolar_lines function from above. You need to do this for both the matrix and the warrior images. A sample output will look as below:
def compute_matching_homographies(e2, F, im2, points1, points2):
'''This function computes the homographies to get the rectified images
input:
e2--> epipole in image 2
F--> the Fundamental matrix (Think about what you should be passing F or F.T!)
im2--> image2
points1 --> corner points in image1
points2--> corresponding corner points in image2
output:
H1--> Homography for image 1
H2--> Homography for image 2
'''
# calculate H2
width = im2.shape[1]
height = im2.shape[0]
T = np.identity(3)
T[0][2] = -1.0 * width / 2
T[1][2] = -1.0 * height / 2
e = T.dot(e2)
e1_prime = e[0]
e2_prime = e[1]
if e1_prime >= 0:
alpha = 1.0
else:
alpha = -1.0
R = np.identity(3)
R[0][0] = alpha * e1_prime / np.sqrt(e1_prime**2 + e2_prime**2)
R[0][1] = alpha * e2_prime / np.sqrt(e1_prime**2 + e2_prime**2)
R[1][0] = - alpha * e2_prime / np.sqrt(e1_prime**2 + e2_prime**2)
R[1][1] = alpha * e1_prime / np.sqrt(e1_prime**2 + e2_prime**2)
f = R.dot(e)[0]
G = np.identity(3)
G[2][0] = - 1.0 / f
H2 = np.linalg.inv(T).dot(G.dot(R.dot(T)))
# calculate H1
e_prime = np.zeros((3, 3))
e_prime[0][1] = -e2[2]
e_prime[0][2] = e2[1]
e_prime[1][0] = e2[2]
e_prime[1][2] = -e2[0]
e_prime[2][0] = -e2[1]
e_prime[2][1] = e2[0]
v = np.array([1, 1, 1])
M = e_prime.dot(F) + np.outer(e2, v)
points1_hat = H2.dot(M.dot(points1.T)).T
points2_hat = H2.dot(points2.T).T
W = points1_hat / points1_hat[:, 2].reshape(-1, 1)
b = (points2_hat / points2_hat[:, 2].reshape(-1, 1))[:, 0]
# least square problem
a1, a2, a3 = np.linalg.lstsq(W, b)[0]
HA = np.identity(3)
HA[0] = np.array([a1, a2, a3])
H1 = HA.dot(H2).dot(M)
return H1, H2
#normalization
def normalize(img):
imgMax = np.amax(img)
imgMin = np.amin(img)
return (img-imgMax)/(imgMax-imgMin)
# convert points from euclidian to homogeneous
def to_homog(points): #here always remember that points is a 3x4 matrix
points2homog = np.vstack((points, [1]*points.shape[1]))
return points2homog
# convert points from homogeneous to euclidian
def from_homog(points_homog):
pointsFhomog = points_homog[:-1,:]/points_homog[-1,:]
return pointsFhomog
def rectify_image(H, source_image):
target_image = np.ones(source_image.shape)
for i in range(target_image.shape[1]):
for j in range(target_image.shape[0]):
ori_points = from_homog(np.matmul(H, to_homog(np.array([[i],[j]]))))
if (0<=ori_points[1]<source_image.shape[0]) and (0<=ori_points[0]<source_image.shape[1]):
target_image[j,i,:] = source_image[int(ori_points[1]),int((ori_points[0])),:]
return target_image
def rectify_points(points,H):
for i in range(points.shape[1]):
points[:, i] = np.dot(H, points[:, i])
# convert the points to cartesian
points[:, i] = points[:, i]/points[:, i][-1]
return points
def image_rectification(im1,im2,points1,points2):
'''this function provides the rectified images along with the new corner points as outputs for a given pair of
images with corner correspondences
input:
im1--> image1
im2--> image2
points1--> corner points in image1
points2--> corner points in image2
outpu:
rectified_im1-->rectified image 1
rectified_im2-->rectified image 2
new_cor1--> new corners in the rectified image 1
new_cor2--> new corners in the rectified image 2
'''
"your code here"
F = fundamental_matrix(points1, points2)
e1, e2 = compute_epipole(F)
H1, H2 = compute_matching_homographies(e2, F.T, im2, points1.T, points2.T)
new_cor1 = rectify_points(points1,H1)
new_cor2 = rectify_points(points2,H2)
rectified_im1 = rectify_image(np.linalg.inv(H1), normalize(im1))
rectified_im2 = rectify_image(np.linalg.inv(H2), normalize(im2))
return rectified_im1,rectified_im2,new_cor1,new_cor2
I1 = imageio.imread("./p4/matrix/matrix0.png")
I2 = imageio.imread("./p4/matrix/matrix1.png")
cor1 = np.load("./p4/matrix/cor1.npy")
cor2 = np.load("./p4/matrix/cor2.npy")
I1, I2, cor1, cor2 = image_rectification(I1, I2, cor1, cor2)
plot_epipolar_lines(I1,I2,cor1,cor2)
I3 = imageio.imread("./p4/warrior/warrior0.png")
I4 = imageio.imread("./p4/warrior/warrior1.png")
cor3 = np.load("./p4/warrior/cor1.npy")
cor4 = np.load("./p4/warrior/cor2.npy")
I1, I2, cor1, cor2 = image_rectification(I3, I4, cor3, cor4)
plot_epipolar_lines(I3,I4,cor3,cor4)
Matching Using epipolar geometry[4 pts]¶
We will now use the epipolar geometry constraint on the rectified images and updated corner points to build a better matching algorithm. First, detect 10 corners in Image1. Then, for each corner, do a linesearch along the corresponding parallel epipolar line in Image2. Evaluate the NCC score for each point along this line and return the best match (or no match if all scores are below the NCCth). R is the radius (size) of the NCC patch in the code below. You do not have to run this in both directions. Show your result as in the naive matching part. Execute this for the warrior and matrix images (Total two outputs images).
def display_correspondence(img1, img2, corrs):
"""Plot matching result on image pair given images and correspondences
Args:
img1: Image 1.
img2: Image 2.
corrs: Corner correspondence
"""
"""
Your code here.
You may refer to the show_matching_result function
"""
show_matching_result(img1, img2, corrs)
def correspondence_matching_epipole(img1, img2, corners1, F, R, NCCth):
"""Find corner correspondence along epipolar line.
Args:
img1: Image 1.
img2: Image 2.
corners1: Detected corners in image 1.
F: Fundamental matrix calculated using given ground truth corner correspondences.
R: NCC matching window radius.
NCCth: NCC matching threshold.
Returns:
Matching result to be used in display_correspondence function
"""
"""
Your code here.
"""
matching = []
corners1 = to_homog(corners1.T)
for i in range(corners1.shape[1]):
a,b,c = np.dot(F.T, corners1[:, i])
max_temp = -10
cmax=()
for x in range(R, img2.shape[1]-R):
y = int((a*x + c)/(-b))
if y >= R and y <= img2.shape[0]-R:
c1 = corners1[0:2, i].astype(np.int32)
c2 = np.array([x, y]).astype(np.int32)
score = ncc_match(img1, img2, c1, c2, R)
if max_temp<score:
max_temp = score
cmax=c2
if max_temp >= NCCth:
matching.append((c1,cmax))
return matching
I1 = imageio.imread("./p4/matrix/matrix0.png")
I2 = imageio.imread("./p4/matrix/matrix1.png")
cor1 = np.load("./p4/matrix/cor1.npy")
cor2 = np.load("./p4/matrix/cor2.npy")
I3 = imageio.imread("./p4/warrior/warrior0.png")
I4 = imageio.imread("./p4/warrior/warrior1.png")
cor3 = np.load("./p4/warrior/cor1.npy")
cor4 = np.load("./p4/warrior/cor2.npy")
# For matrix
rectified_im1,rectified_im2,new_cor1,new_cor2 = image_rectification(I1,I2,cor1,cor2)
F_new = fundamental_matrix(new_cor1, new_cor2)
nCorners = 10
# Choose your threshold
NCCth = 0.7
#decide the NCC matching window radius
R = 25
# detect corners using corner detector here, store in corners1
corners1 = corner_detect(rgb2gray(rectified_im1), nCorners, smoothSTD, windowSize)
corrs = correspondence_matching_epipole(rgb2gray(rectified_im1), rgb2gray(rectified_im2), corners1, F_new, R, NCCth)
display_correspondence(rectified_im1, rectified_im2, corrs)
# For warrior
rectified_im3,rectified_im4,new_cor3,new_cor4 = image_rectification(I3,I4,cor3,cor4)
F_new2=fundamental_matrix(new_cor3, new_cor4)
nCorners = 10
# Choose your threshold
NCCth = 0.7
#decide the NCC matching window radius
R = 25
# You may wish to change your NCCth and R for warrior here.
corners2 = corner_detect(rgb2gray(rectified_im3), nCorners, smoothSTD, windowSize)
corrs = correspondence_matching_epipole(rgb2gray(rectified_im3), rgb2gray(rectified_im4), corners2, F_new2, R, NCCth)
display_correspondence(rectified_im3, rectified_im4, corrs)
